using Distributions,Printf,Random,Statistics

#const α    = 2.0
#const β    = 5.0
const α    = 2.5
const β    = 5.0

const grid = 0.0000001
const M    = 10000000
const ρ  = [-0.999,-0.99,-0.9,-0.5,0.0,0.5,0.9,0.99,0.999] 
results11  = Array{Float64,2}(undef,length(ρ),11) 
results12  = Array{Float64,2}(undef,length(ρ),11) 
results21  = Array{Float64,2}(undef,length(ρ),11) 
results22  = Array{Float64,2}(undef,length(ρ),11) 
Random.seed!(2024)

# Sign restrictions on demand and supply

for i=1:length(ρ)
    φ = asin(ρ[i])
    if φ<0.0
        μ1 = pi/2.0-φ
        μ2 = pi
    else
        μ1 = pi/2.0
        μ2 = pi-φ 
    end
    results11[i,1] = minimum(map(μ->cos(μ),collect(μ1:grid:μ2)))
    results11[i,2] = maximum(map(μ->cos(μ),collect(μ1:grid:μ2)))
    results12[i,1] = minimum(map(μ->sin(μ),collect(μ1:grid:μ2)))
    results12[i,2] = maximum(map(μ->sin(μ),collect(μ1:grid:μ2)))
    results21[i,1] = minimum(map(μ->sin(μ+φ),collect(μ1:grid:μ2)))
    results21[i,2] = maximum(map(μ->sin(μ+φ),collect(μ1:grid:μ2)))
    results22[i,1] = minimum(map(μ->-cos(μ+φ),collect(μ1:grid:μ2)))
    results22[i,2] = maximum(map(μ->-cos(μ+φ),collect(μ1:grid:μ2)))    
    # prior 1
    μ = 2.0*pi*rand(M)
    μ = μ[(μ1.<=μ).&(μ.<=μ2)]
    N = length(μ) 
    display(N)
    A = [cos.(μ) sin.(μ) sin.(φ*ones(N)+μ) -cos.(φ*ones(N)+μ)]
    A = sort(A,dims=1) 
    results11[i,3] = mean(A[:,1])
    results12[i,3] = mean(A[:,2])
    results21[i,3] = mean(A[:,3])
    results22[i,3] = mean(A[:,4])
    results11[i,4] = A[Int64(floor(0.05*N)),1] 
    results12[i,4] = A[Int64(floor(0.05*N)),2] 
    results21[i,4] = A[Int64(floor(0.05*N)),3] 
    results22[i,4] = A[Int64(floor(0.05*N)),4] 
    results11[i,5] = A[Int64(floor(0.95*N)),1] 
    results12[i,5] = A[Int64(floor(0.95*N)),2] 
    results21[i,5] = A[Int64(floor(0.95*N)),3] 
    results22[i,5] = A[Int64(floor(0.95*N)),4] 

    # prior 2
    μ = 2.0*pi*rand(Beta(α,β),M)
    μ = μ[(μ1.<=μ).&(μ.<=μ2)]
    N = length(μ) 
    display(N)
    A = [cos.(μ) sin.(μ) sin.(φ*ones(N)+μ) -cos.(φ*ones(N)+μ)]
    A = sort(A,dims=1) 
    results11[i,6] = mean(A[:,1])
    results12[i,6] = mean(A[:,2])
    results21[i,6] = mean(A[:,3])
    results22[i,6] = mean(A[:,4])
    results11[i,7] = A[Int64(floor(0.05*N)),1] 
    results12[i,7] = A[Int64(floor(0.05*N)),2] 
    results21[i,7] = A[Int64(floor(0.05*N)),3] 
    results22[i,7] = A[Int64(floor(0.05*N)),4] 
    results11[i,8] = A[Int64(floor(0.95*N)),1] 
    results12[i,8] = A[Int64(floor(0.95*N)),2] 
    results21[i,8] = A[Int64(floor(0.95*N)),3] 
    results22[i,8] = A[Int64(floor(0.95*N)),4] 

    # prior 3
    μ = 2.0*pi*rand(Beta(β,α),M)
    μ = μ[(μ1.<=μ).&(μ.<=μ2)]
    N = length(μ) 
    display(N)
    A = [cos.(μ) sin.(μ) sin.(φ*ones(N)+μ) -cos.(φ*ones(N)+μ)]
    A = sort(A,dims=1) 
    results11[i,9] = mean(A[:,1])
    results12[i,9] = mean(A[:,2])
    results21[i,9] = mean(A[:,3])
    results22[i,9] = mean(A[:,4])
    results11[i,10] = A[Int64(floor(0.05*N)),1] 
    results12[i,10] = A[Int64(floor(0.05*N)),2] 
    results21[i,10] = A[Int64(floor(0.05*N)),3] 
    results22[i,10] = A[Int64(floor(0.05*N)),4] 
    results11[i,11] = A[Int64(floor(0.95*N)),1] 
    results12[i,11] = A[Int64(floor(0.95*N)),2] 
    results21[i,11] = A[Int64(floor(0.95*N)),3] 
    results22[i,11] = A[Int64(floor(0.95*N)),4] 
end
io = open("uhlig_both.tex","w")
@printf(io,"\\documentclass[12pt]{article}\n")
@printf(io,"\\usepackage{rotating}\n")
@printf(io,"\\begin{document}\n")
@printf(io,"\\footnotesize\n")
@printf(io,"\\begin{sidewaystable}\n")
@printf(io,"Price response to supply shock\n")
@printf(io,"\\vspace*{1em}\n\n")
@printf(io,"\\begin{tabular}{l|ll|lll|lll|lll}\\hline\n")
@printf(io," & & & \\multicolumn{9}{|c}{posterior for \$\\theta\$} \\\\ \\cline{4-12} \n")
@printf(io," & & & \\multicolumn{3}{|c|}{using prior 1} &  \\multicolumn{3}{|c|}{using prior 2} & \\multicolumn{3}{|c}{using prior 3} \\\\ \\cline{4-12}\n")
@printf(io," \$\\rho\$ & \\multicolumn{2}{|c|}{identified set} & mean & \\multicolumn{2}{c|}{90\\%% CI} & mean & \\multicolumn{2}{c|}{90\\%% CI} & mean & \\multicolumn{2}{c}{90\\%% CI} \\\\ \\hline \n")
for i=1:length(ρ)
    @printf(io,"%5.3f ",ρ[i])
    @printf(io,"& [%6.3f, ",results11[i,1])
    @printf(io,"& %6.3f] ",results11[i,2])
    for j=1:3
        @printf(io,"& %6.3f ",results11[i,3*j])
        @printf(io,"& [%6.3f, ",results11[i,3*j+1])
        @printf(io,"& %6.3f] ",results11[i,3*j+2])
    end
    @printf(io,"\\\\ \\hline \n") 
end
@printf(io,"\\end{tabular}\n\n")
@printf(io,"\\vspace*{2em}\n")
@printf(io,"Price response to demand shock\n")
@printf(io,"\\vspace*{1em}\n\n")
@printf(io,"\\begin{tabular}{l|ll|lll|lll|lll}\\hline\n")
@printf(io," & & & \\multicolumn{9}{|c}{posterior for \$\\theta\$} \\\\ \\cline{4-12} \n")
@printf(io," & & & \\multicolumn{3}{|c|}{using prior 1} &  \\multicolumn{3}{|c|}{using prior 2} & \\multicolumn{3}{|c}{using prior 3} \\\\ \\cline{4-12}\n")
@printf(io," \$\\rho\$ & \\multicolumn{2}{|c|}{identified set} & mean & \\multicolumn{2}{c|}{90\\%% CI} & mean & \\multicolumn{2}{c|}{90\\%% CI} & mean & \\multicolumn{2}{c}{90\\%% CI} \\\\ \\hline \n")
for i=1:length(ρ)
    @printf(io,"%5.3f ",ρ[i])
    @printf(io,"& [%6.3f, ",results12[i,1])
    @printf(io,"& %6.3f] ",results12[i,2])
    for j=1:3
        @printf(io,"& %6.3f ",results12[i,3*j])
        @printf(io,"& [%6.3f, ",results12[i,3*j+1])
        @printf(io,"& %6.3f] ",results12[i,3*j+2])
    end
    @printf(io,"\\\\ \\hline \n") 
end
@printf(io,"\\end{tabular}\n")
@printf(io,"\\end{sidewaystable}\n")
@printf(io,"\\clearpage\n")
@printf(io,"\\begin{sidewaystable}\n")
@printf(io,"Quantity response to supply shock\n")
@printf(io,"\\vspace*{1em}\n\n")
@printf(io,"\\begin{tabular}{l|ll|lll|lll|lll}\\hline\n")
@printf(io," & & & \\multicolumn{9}{|c}{posterior for \$\\theta\$} \\\\ \\cline{4-12} \n")
@printf(io," & & & \\multicolumn{3}{|c|}{using prior 1} &  \\multicolumn{3}{|c|}{using prior 2} & \\multicolumn{3}{|c}{using prior 3} \\\\ \\cline{4-12}\n")
@printf(io," \$\\rho\$ & \\multicolumn{2}{|c|}{identified set} & mean & \\multicolumn{2}{c|}{90\\%% CI} & mean & \\multicolumn{2}{c|}{90\\%% CI} & mean & \\multicolumn{2}{c}{90\\%% CI} \\\\ \\hline \n")
for i=1:length(ρ)
    @printf(io,"%5.3f ",ρ[i])
    @printf(io,"& [%6.3f, ",results21[i,1])
    @printf(io,"& %6.3f] ",results21[i,2])
    for j=1:3
        @printf(io,"& %6.3f ",results21[i,3*j])
        @printf(io,"& [%6.3f, ",results21[i,3*j+1])
        @printf(io,"& %6.3f] ",results21[i,3*j+2])
    end
    @printf(io,"\\\\ \\hline \n") 
end
@printf(io,"\\end{tabular}\n\n")
@printf(io,"\\vspace*{2em}\n")
@printf(io,"Quantity response to demand shock\n")
@printf(io,"\\vspace*{1em}\n\n")
@printf(io,"\\begin{tabular}{l|ll|lll|lll|lll}\\hline\n")
@printf(io," & & & \\multicolumn{9}{|c}{posterior for \$\\theta\$} \\\\ \\cline{4-12} \n")
@printf(io," & & & \\multicolumn{3}{|c|}{using prior 1} &  \\multicolumn{3}{|c|}{using prior 2} & \\multicolumn{3}{|c}{using prior 3} \\\\ \\cline{4-12}\n")
@printf(io," \$\\rho\$ & \\multicolumn{2}{|c|}{identified set} & mean & \\multicolumn{2}{c|}{90\\%% CI} & mean & \\multicolumn{2}{c|}{90\\%% CI} & mean & \\multicolumn{2}{c}{90\\%% CI} \\\\ \\hline \n")
for i=1:length(ρ)
    @printf(io,"%5.3f ",ρ[i])
    @printf(io,"& [%6.3f, ",results22[i,1])
    @printf(io,"& %6.3f] ",results22[i,2])
    for j=1:3
        @printf(io,"& %6.3f ",results22[i,3*j])
        @printf(io,"& [%6.3f, ",results22[i,3*j+1])
        @printf(io,"& %6.3f] ",results22[i,3*j+2])
    end
    @printf(io,"\\\\ \\hline \n") 
end
@printf(io,"\\end{tabular}\n\n")
@printf(io,"\\end{sidewaystable}\n")
@printf(io,"\\end{document}\n")
close(io)

